Spikes in the Mixmaster regime of G2 cosmologies 
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I. INTRODUCTION 

Belinskii, Khalatnikov and Lifshitz (BKL) con- 
jectured that, according to general relativity, the ap- 
proach to the generic spacelike singularity is vacuum 
dominated (assuming p < p), local, and oscillatory (la- 
beled 'Mixmaster'). Here local means that the contri- 
bution of terms in the evolution equations with spatial 
derivatives becomes negligible. The Mixmaster dynam- 
ics consists of Kasner epochs bridged by transitions rep- 
resented by the vacuum Bianchi type II solution. Nu- 
merical studies of the asymptotics of the Gowdy mod- 
els, which represent the simplest inhomogeneous vacuum 
spacetimes, reveal that on approach to the singularity, 
spiky structures form 0, S [B| • These spikes become ever 
narrower as the singularity is approached. At first, the 
presence of such spikes might seem inconsistent with the 
'local' part of the BKL conjecture since the spatial deriva- 
tive of such a spike grows without bound as the singu- 
larity is approached. Remarkably, such spiky behavior in 
Gowdy spacetimes actually is consistent with BKL local- 
ity. The reason for this is that in the evolution equations 
for Gowdy spacetimes, the spatial derivatives are mul- 
tiplied by a quantity that goes to zero even faster than 
the spatial derivatives go to infinity. This has been ver- 
ified in detailed numerical simulations, as well as by the 
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discovery of closed form Gowdy solutions with the spike 
property 

Nonetheless, the Gowdy models are a very special class 
of spacetimes, so it remains to be seen whether this prop- 
erty of BKL locality and persistent spikes holds in more 
general spacetimes. To that end, we will examine the 
properties of spikes in G2 models which are a slightly 
more general class that includes the Gowdy spacetimes. 

Studies of G2 and more general models have pro- 
duced numerical evidence that the BKL conjecture gener- 
ally holds except poss ibly at isolated points where spiky 
structures form [1 HO, El [3 • Here, BKL locality IS vio- 
lated due to large spatial gradients. However, the ability 
to draw conclusions about spikes from such simulations is 
severely limited due to the enormous numerical resources 
needed to resolve the narrowing spikes. In this paper, we 
will use a different numerical method that does have ad- 
equate resolution to provide reliable conclusions about 
spike behavior in G2 spacetimes. We present numerical 
evidence in support of the following conjectures: that 
recurring "spike transitions" are a general type of oscil- 
lation as the singularity is approached, and that higher- 
order spike transitions split into separate first-order spike 
transitions. Section 2 presents the equations for the evo- 
lution of G2 spacetimes. Our numerical method is pre- 
sented in section 3, results in section 4, and conclusions in 
section 5. Appendix [K\ presents the procedure to match 
a numerical solution with an explicit spike solution. Ap- 
pendix |B] gives the formula for the BKL parameter u in 
term of the parameter w. Appendix[C] gives the formulae 
for the Weyl scalar invariants. 



2 



II. G2 SPACETIMES 

The metric of the general G2 class takes the form [13, 
eq (7)] 

+ e^-" [dy + g dz + (Gi + QG2)dx 

+ (A/i + gAf2)(-e-Mr)]2 

+ e-^-"[dz + G2dx + A/2(-e-"dr)]2. (1) 

Here all metric quantities depend only on the time coor- 
dinate r and spatial coordinate x, thus there is symmetry 
in two spatial directions. The singularity is approached 
as T ^ c». The choice of gauge used here is the same as 

in 

Our choice of variables are the /3-normalized variables 
[13I [T3 | in the orthonormal frame formalism J^] , related 
to the metric components as follows: 

1 

- AA=--, (2) 



f3 = 2e- 
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e-4— 



V3 ' 



(3) 
(4) 
(5) 



73 ' V3 ' 

where if is a constant, and the t and a; subscripts denote 
partial differentiation. 

The evolution equations for the /3-normalized variables 



drEi^ ^ -\(2 ~ i-El)E^^ (6) 
a^E- = -i[-3I]2S_ + 2\/3(S2 - Ar2) _ ^5^2 

-E^^d^N^] (7) 
5riVx = -i[(2 - 3Sl)iVx - (8) 
a^Ex = -i[(-3I]^ - 2V3E_)Ex - 2V3A^xA^- 

+ E^^d,,N^] (9) 
drN_ = - i[(2 - 3E^ + 2V3S_)iV_ + 2V3SxiVx 

+ E^^d^^y.] (10) 

a^Ea = -5[-3E^-3E+ + V3I]_]E2, (11) 
where 

1 



-(1 - - E^ - E^ - A^l - Nl). 



(12) 



There is one constraint equation: 

Ei^d^^2 = (3A^-Ex - 3A^xE_ - V3iVx)E2. (13) 

For state-space presentations, we will use the Hubble- 
normalized variables @: 

(E+,E_,Ex,E2,iV_,iVx)^ 
1 



1-E. 



(E+,E_,Ex,E2,7V-,7Vx). (14) 



See |12[ for the evolution equations for Hubble- 
normalized variables, and the derivation of the evolution 
equations. 

The Gowdy spacetimes are that class of G2 spacetimes 
for which E2 = 0. Note that it then follows from equation 
^ that Ei^ — exp(— r). An interesting class of solutions 
of the Gowdy equations are the exact spike solutions of 



E_ = (1 + [w tanh(,«T) - 1]) (15) 



A^x 
Ex 



2/ w 
/2 + I V3 

P -I w 

2/ 1 



sech (wt) 
sech {wt) 
{1 — w tanh(wr)) 



where w is a constant and the quantity / is given by 



(16) 
(17) 
(18) 



/ = we^ sech {wt)x 



(19) 



For |w| < 1 this solution describes a spike because f — 
at a; = but / becomes large as r — 00 for all a: ^ 0. 
Nonetheless, BKL locality is preserved because in the 
equations of motion all spatial derivatives are multiplied 
by Ei^ and Ei^dxf = w sech (wt) which goes to zero as 
r 00. 

Note however that this conclusion depends on the fact 
that Ei^ — exp(— r) which in turn depends on the fact 
that in equation ^ we could set E2 to zero, something 
that we can only do in Gowdy spacetimes, not the more 
general G2 spacetimes. The dynamics in a G2 spacetime 
consists of eras where E2 is very small (and which can 
thus be well described by the dynamics of Gowdy space- 
times) punctuated by short "frame bounces" where E2 
rapidly grows and then rapidly shrinks to become again 
negligible. During a frame bounce Ei^ shrinks more 
slowly than exp(— r) and thus it is not clear whether spa- 
tial derivatives continue to remain negligible. To resolve 
this issue, we will need to perform numerical simulations 
of the dynamics of G2 spacetimes. Furthermore, those 
simulations will need to have enough resolution to accu- 
rately model the rapidly shrinking spikes. 



III. NUMERICAL METHODS 

One numerical method for resolving small scale struc- 
ture is adaptive mesh refinement (AMR). However, if one 
knows beforehand the location of the structure, one need 
not use AMR and can instead use a coordinate system 
adapted to the structure that one wants to study. In 
particular, here we are studying spikes that shrink ex- 
ponentially with time, so we choose a coordinate system 
that docs the same. 

We introduce new coordinates {T,X) to zoom in on 
the worldline x = a;,norTi. 



T 



X 



{x a^zoom). 



(20) 
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FIG. 1: The spacetime diagram showing the x — const world- 
hnes (vertical), the spatial hypersurfaces t = T — const 
(horizontal) with t = T cx> at the singularity, the par- 
ticle horizon (45° lines) of the observer in the center, and the 
X = const lines which are timelike inside the particle horizon 
but spacelike outside (dashed lines). 



where the constant A controls the rate of focus. See 
Figure [T] for the qualitative spacetime diagram. The dif- 
ferential operators expressed in the new coordinates are 



dr = dT + AXdx, d^^e^^dx 



The equations in the new coordinates are 



(21) 



St El ^ 
9x^2 = — 



-AXdxEi^ - [1 



-AXdx^- + ^e^^Ei^dxNx 



\/3 v2 
— ^2 



(22) 



(23) 



- \/3(E^ - 

^AXdxNx + \e^^Ei^dx^- - [1 - fS^lA^x 

(24) 

-AXdx^x - ^e^'^Ei^dxN^ + fS^Ex 

(25) 

(26) 
(27) 



-AXdxN_-^e^^Ei^dx^x-[l 
AXdx^2 + [fSa + 1^+ - ^S- 



and constraint 



e'^^Ei^dx^2 - (37V_Ex - 3iVxS- - V3iVx)S2. (28) 



We will end the numerical grid at a fixed coordinate 
value X — Xq . Ordinarily, that would call for a boundary 
condition at Xq, but we will use the method of excision. 
Usually one thinks of excision as applying to simulations 
of black holes; however excision can be applied to any hy- 
perbolic equations where the outer boundary is chosen so 
that all modes are outgoing. In that case one simply im- 
plements the equations of motion at the outer boundary, 
no boundary condition is needed (or even allowed). 



The following combinations of the equations of motion 
9t(S_ +iVx) = -{AX - ^e"^^ Ei^)dxiy'- + Nx) + ■ ■ ■ 



(AX + ^e^^Ei')dx{^--Nx 



aT(S- -iVx) 
(9t(Sx +iV-) = -{AX + ^c'^^Ei^)dx{T:x +iV_ 



-{AX-^c^^Ei')dx{^x-N^ 



(29) 
^30) 

^31) 

1 

(32) 



clearly shows that (S- + N^) and (Sx — flow away 
from X = jje-^^Ei^ (for A > 0), while (S_ - A^x) and 
(Ex + iV_) flow away from X = -^e'^'^Ei'^. This puts 
X — ±j^e^^Ei^ as the points beyond which the flow 
is entirely outward. Thus, as long as Xq is chosen large 
enough and as long as e^^Ei^ does not grow too large 
during the simulation, the surface X = Xq will be a good 
excision boundary. 

In addition to choosing Xq, we should also choose A so 
that X = Xq remains a good excision boundary through- 
out the simulation. A — I is the natural choice, which 
fixes the particle horizon of the exact spike solution as 
a vertical line in the spacetime diagram with respect to 
{X,T) coordinates. In this paper we shall choose A = 1. 
Choosing another value for A is a trial and error process, 
but one is able estimate Ei^ after one or two numerical 
runs, with the heuristics below. 

We shall define phenomenologically that a Gowdy era 
as the time period during which E2 is small. We take this 
opportunity to correct that the Kasner eras mentioned in 

are in fact Gowdy eras. The two are not equivalent, as 
there can be two or three Kasner eras within one Gowdy 
era. During a Gowdy era, Ei^ approximately equals e~'^, 
but between Gowdy eras (namely during the E2 transi- 
tion) Ei^ shrinks more slowly. In order to offset this 
behavior between Gowdy eras, one should choose a small 
enough A < 1 so that e^^Ei^ decays during the Gowdy 
era. But if A is too small, spikes will be inadequately 
resolved. A reasonable range is 0.8 < A < 1. 

Another way to make X = Xq a good excision bound- 
ary is to choose a larger Xq to leave more room for the 
growth of e^^Ei^. The CFL condition, however, requires 
that the numerical time step AT satisfies 



AT < AXq 



^AT 1 



El'] AX, 



(33) 



where AX is the numerical grid size. For example, dou- 
bling Xq would cut AT by a half, so one is bound by nu- 
merical resources to choose a large enough Xq for the sim- 
ulation without being too wasteful. A reasonable range 
is 10 < ATq < 40. 

Our numerical simulations use a uniform spatial grid. 
The equations are evolved using the classical fourth-order 
Runge-Kutta method, with fourth-order accurate spatial 
derivatives. That is for any quantity F we approximate 



4 



dxF on grid point i by 



4F,+i 
3" 



F,. 



F,, 



2AX 



l F,+2 

3 4AX 



i-2 



(34) 



On the last gridpoint, the excision boundary, we evaluate 
the spatial derivative using one sided differences. That 
is we approximate dxF at the final gridpoint N by 



25Fjv - 48j^jv„i + 36fAr^2 - 16f jv-3 + 3f jv-4 
12AX 

and at the second last grid point iV — 1 by 

3f AT + IOFat-I - 18FAf,2 + 6FAf-3 - -FjV-4 



(35) 



(36) 



For spikes, non-symmetric data would be problematic for 
implementing the local perspective as the spike worldline 
is not stationary in this case. Therefore we shall choose 
symmetric initial data (around X = Q) and simulate only 
X e [0,Xo], with enforcement of the symmetry at the 
left boundary X = 0. For comparison, we also simulate 
along non-spike worldlines, in which case the data are not 
symmetric and the left boundary at —X^ is an excision 
boundary. 

We choose the first gridpoint to be either an excision 
boundary &t X = —Xq or a point of symmetry dX X — Q. 
If it is an excision boundary, then dxF is approximated 
by the one sided differences 



-Fi + 48F2 - 36F3 + I6F4 - 3i^5 



12AX 



at the first grid point, and 



-3i^i - IOF2 -f 18i^3 - 6F4 + F5 
12AX 



(37) 



(38) 



at the second. However, if first grid point is a point of 
symmetry then we choose all quantities to be either even 
or odd there. For even functions, dxF = a,t the first 
grid point, and 



4 — Fi 1 F4 — F2 



3 2AX 3 AAX 



(39) 



and the second, while for odd functions we approximate 
dxFhy 



SF2 — i^3 

6AX 

at the first grid point, and 



4 F3 - j^i 
3 2AX 



IF4 

3 4AX 



F2 



(40) 



(41) 



at the second. 

The standard double precision real variables (with 16 
digits of significance) are normally used in the numeri- 
cal code. When necessary, quad precision real variables 
(with 32 digits of significance) are used to lower the nu- 
merical roundoff errors by 10^^ folds, thereby preventing 



it from prematurely swamping small values. The vari- 
ables iV_, Sx and A^x take small values during Kasner 
epochs, and can be swamped by the roundoff error in 
the spatial derivative term of another variable with a 
larger value. Usually this happens to iVx first, when the 
term ^e^^ Ei^dx^- in equation ([M)) becomes 10^^ times 
smaller (if double precision is used) than the value of 
The usage of quad precision real variables increases the 
runtime by 4 to 8 folds. 

No numerical dissipation is used, as it is unnecessary. 

To verify that numerical solutions converge with fourth 
order accuracy, we compare the constraint (j27p in numer- 
ical runs with different resolutions (different number of 
grid points). We observe that doubling the resolution re- 
duces the constraint by a factor of 16 when adequate nu- 
merical resolution is used. We also compare the numer- 
ical solutions with a matching exact spike solution. The 
procedure for matching is described in Appendix El The 
formula for the BKL parameter u for the Kasner epochs 
between transitions are given in Appendix |B] The Weyl 
scalar invariants are used to measure the difference be- 
tween numerical and exact solutions. Appendix [C] gives 
their formulae. 

In this paper, we shall focus on obtaining numerically 
accurate results, which require much higher numerical 
resolution than qualitative numerical results do. This re- 
quirement also places severe limit on how far into the 
asymptotic regime one can simulate, because the numer- 
ical error must not be larger than the distance from the 
solution to the nearest Kasner point in the state space, 
and this in turn require high numerical resolution. When 
a numerical simulation takes up to months to run in or- 
der to meet the accuracy, it becomes impractical. Despite 
this difficulty, we want to provide more than just qualita- 
tive numerical results, because numerically accurate re- 
sults can provide evidence supporting convergence to the 
exact spike solution, while qualitative numerical results 
cannot. In presentation, we shall round the numbers to 
4 decimals, even though the accuracy is higher. 

Qualitative numerical results are still valuable in pro- 
viding evidence supporting the general behavior of the 
solution. Compared with other aspects of the solution, 
the timing of a transition is most sensitive to numerical 
inaccuracy. At lower resolutions, the timing of a tran- 
sition differs greatly while other aspects of the solution 
remain robust. 



IV. RESULTS 

We clarify a few terms we use below. A (true) spike 
point is where a (true) spike can occur (the spike may be 
active or smoothed out). In our variables, a spike point 
is where 



7V_ = 0. 



(42) 



We shall hold the spike point fixed (at a; = 0), so that 
we can easily locate it and zoom in on it. To do so, we 
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require iV_ and to be odd functions around the spike 
point, and S]_, Sx and S2 to be even functions. 

A false spike point is where Sx = 0. To hold it fixed, 
we require Sx and A^x to be odd functions around the 
false spike point, and S-, S2 and N- to be even func- 
tions. 

We will present three sets of numerical results. The 
first set chooses a perturbed spike solution as the initial 
condition, and shows two recurrences of the spike solu- 
tion, within the same Gowdy era. The purpose is to show 
that spike recurs within the same Gowdy era. The sec- 
ond set chooses a generic initial condition, and shows two 
occurrences of the spike solution, one in each Gowdy era. 
The purpose is to show that spike recurs over different 
Gowdy eras. The third set consists of two simulations, 
with a perturbed second and third order spike solution 
as the initial condition, respectively. The purpose is to 
show that second and third order spikes break up into 
separate first order spikes. All three sets demonstrate 
the attractor nature of the first order spike. 

The reader will notice that different numerical resolu- 
tions are used for different sets. The length of simulation 
also differs. Both numerical resolutions and length are 
not arbitrarily chosen, but are dictated by the cost of 
computation to maintain accuracy. For example, in the 
first set, we do not try to extend the simulation to show 
the third spike recurrence over different Gowdy eras, as 
it would be too costly. 



A. Perturbed spike 

The format of the initial data is a perturbed spike so- 
lution at r = 0: 



El' = 2, 
11 1 



Sx = 



{wxf + 1^/3 \/3 
(wx)'^ — 1 w 



S2 = 10 



A^. 



-5 



{wxy + 1 Vs' 



(43) 

2wx w 

{wx)^ + 1 Vl' 
(44) 

2wx 1 

(45) 



(wx)^ + 1 VS' 



where x = X + x^oomj with e = 0.001, w = 9.5, A — 
1, a^zoom — 0. Here, small perturbations are applied to 
the variables S2 and E_ . Compare with the exact spike 
solution at T = [1, eq (36)]. The value w = 9.5 allows 
for two spike recurrences at roughly w = 5.5 and w = 1.5 
before the next Gowdy era. 

Because the data chosen are symmetric about X ~ 
0, only A > needs to be simulated. A resolution of 
100001 grid points on the A'-interval [0,2] is sufficient for 
convergence during the time interval T S [0, 16]. Double 
precision is used. Beyond T = 16 a higher numerical 
precision is needed to maintain accuracy. To compare 
the orbit along a different worldline, we also use the same 
initial data but with Xzoom = 1 Sind 200001 grid points 
on the A'-interval [—2,2]. 



The simulation shows a perturbed spike solution re- 
curs twice over the same Gowdy era. For each of the 
two recurrences, the numerical solution is matched with 
an exact spike solution. The difference between the nu- 
merical and exact spike solution is computed in the four 
Weyl scalar invariants, and is observed to be smaller in 
the second recurrence than in the first (see Figures [2] and 
131). This suggests that the closer to the singularity, the 
closer the numerical solution gets to an exact spike solu- 
tion. This supports the conjecture that the exact spike 
solutions are attractors. 

Figure d] shows the orbit along the spike point x = 
(from the first simulation) and the worldline x — 1 (from 
the second simulation) projected onto the (E:^,E:^^) 
plane in the state space of Hubble-normalized variables. 
It shows the orbits follow the expected paths as predicted 
from the spike solution and the Ex and A^_ transition 
sets (see Figures 5 and 6 in Q). This subsection is simi- 
lar to the work in , which was done in the Gowdy class 
(E2 = 0), in the sense that the simulations here focus on 
what happens within one Gowdy era. The approximate 
values for the w parameter in [8| and corresponding u 
parameter when near a Kasner point for the x = or- 
bit in Figure [4] are given below (rounded to 4 decimal 
points). Linking the Kasner epochs are alternating frame 
and spike transitions. 



w « -7.8330 7.5401 -3.7592 ''-^^ 3.8264 



spike 



0.1674 (46) 

M-2 



u « 3.4165 3.2700 ^ 1.3796 1.4132 
^ 1.4021 



(47) 



Note that Kasner epochs linked by frame transitions are 
not distinct physically. One also observes that the num- 
bers above do not follow the maps very closely, suggesting 
that the solution is not yet very close to the generalized 
Mixmaster attractor. The difference is due to perturba- 
tion present in the initial data. Over time, the difference 
gradually decreases. Also note that the map u — > m — 2 
has an adjustment algorithm when the new value is less 
than 1, namely 



u — 2 if u > 3 

if 2 < u < 3 
if I < u < 2 



1 

u-2 

1 

~r — 

i-i ■ 



(48) 



B. Generic initial condition 

Having seen that perturbed spike initial data lead to 
recurring spikes within the same Gowdy era, we now go 
further and ask whether a generic initial data also lead to 
recurring spikes, and whether the recurrence continues in 
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Numerical solution 



Exact solution 



Difference 



CC/(3hY 






100 

CCC/{3H^f g 



-100 



100 

CC*C/(3H^)^ g 

-100 






100 


-100 




FIG. 2: The Weyl scalar invariants (normalized by Hubble scalar) for the numerical solution, the matching exact spike solution 
(with w — 5.8644, to = 2.2205), and their difference, during the first spike recurrence over the time interval [0.5, 3.5]. 



the next Gowdy era. To answer these questions, we shall 
start with a generic initial data, and evolve the solution 
through to the next Gowdy era. 

We give an example of a generic initial condition for a 
true spike below. 

S-=ai+a2X^, Nx = a^x, (50) 



Sx = as + a4:X , 



N- — a^x, 



(51) 



with 



/ Cx/Ei'dx — -(— Sflifls — VSi 
Jo 4 



0-5 + 30306)2;^ 



1 



o 



where x = X + Xy, 



For example, we choose 
ai 3.25, aa = 0.002, ag = 0.3, 
04 = -0.001, as = 0.04, as = -0.05, 
(^2)0 = 0.2, 



(52) 

(53) 
(54) 
(55) 



with A = 1, Xzoom = 0, 6401 grid points over the X- 
interval [0, 10], and time interval [0,40]. Quadruple pre- 
cision is used. For comparison, another simulation with 
a^zoom = 1, 12801 grid points over the X-interval [—10, 10] 
is used. Beyond T = 40, the solution gets too close to a 
Kasner point, and a higher numerical resolution is needed 
to maintain accuracy. 



Two recurrences of spike are observed, one in the same 
Gowdy era, and the other in the next (after a S2 tran- 
sition). Figure [5] shows the orbits along a; = 0, 1 passing 
close to various identical points during two Gowdy eras. 
A difference in position of the final points is observed, 
and is attributed to the lag between the two worldlines 
that becomes more pronounced over time. The approx- 
imate values for the w and u parameters when near a 
Kasner point (except the initial point) for the x — or- 
bit in Figure are given below (rounded to 4 decimal 
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FIG. 4: The orbit along the spike point x — and a nearby non-spike point x — 1 projected onto the (E+ , E^f) plane for the 
perturbed spike simulation. The initial (at T = 0) and final (at T = 16) points are marked by the letters i and f respectively. 



points). 



w « 5.7070 -L7114 ''-^^ T7114 6.6228 

-2.6228 (56) 

1 

u w 2.3535 ^ 2.8114 2.8114 2.8114 

1 

^ 1.2324 (57) 



Observe that the numbers here follow the map much 
more closely in the later stage. 

The numerical solution is matched with an exact spike 
solution and the Weyl scalars are plotted in Figures [5] 
and [71 As in the previous subsection, matching improves 
with time (towards the singularity). The remarkable im- 
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FIG. 5: The orbit along the spike point x = and a nearby non-spike point x — 1 projected onto the (S+ , E^) plane for the 
generic initial data simulation. The initial (at T = 0) and final (at T — 40) points are marked by the letters i and f respectively. 





FIG. 6: The Weyl scalar invariants (normalized by Hubble scalar) for the numerical solution, the matching exact spike solution 
(with w — —3.7114, To = 2.6461), and their difference, during the first spike recurrence over the time interval [0, 5]. 



provement from the first to the second recurrence also 
suggest an exponential rate of convergence to the exact 
spike solution. This provides a very strong evidence that 
the spike solution is an attractor not only for perturbed 
spike initial data, but also for generic ones. The expo- 
nential rate of convergence is also a curse for accurate 



numerical simulations, as the need for numerical resolu- 
tion also increases exponentially with time. 
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Numerical solution 



Exact solution 



Difference times 10 



11 



CC/(3hV 



C*C/(3hY 



240 

CCC/(3hY 




240 

CC*C/(3hY q 



38-5 
T 36 X 

FIG. 7: The Weyl scalar invariants (normalized by Hubble scalar) for the numerical solution, the matching exact spike solution 
(with w — —4.6228, tq = 38.1923), and their difference, during the second spike recurrence over the time interval [35,40]. 



C. Perturbed higher order spikes 

Having seen that the spike solution is an attractor, we 
now investigate whether higher order spikes 0, Section 
5.5] are also attractors. In this subsection we shall use 
perturbed second and third order spike solutions as initial 
data, and see whether they recur as first, second or third 
order spikes. 

The initial data for a perturbed second order spike 
solution at r = is given recursively in terms of the first 
order spike solution: 



Ei^ = 2, 
1 

"71' 



(58) 

TVx =sAf-i-ciVxi, 
(59) 



Af- = -sS_i +cSxi, (60) 



spike solution in c and s are given by 



h 



1 



-2w(l - 2x^) + 



Q2 



(61) 
(62) 



where x = X + x^oom- 
evaluating the constraint 



2go_ 

E2 is specified by numerically 



S2 = (£2)06^°^^^/^^'''". (63) 
Wc perform a simulation with the parameters 



e = 10"^ w = 9.5, 

(£2)0 - lo-^ 



?2 = 0, 



(64) 
(65) 



where (E_, A^x , Sx , iV_)i are the perturbed first order Te [0,12]. 



0, X e [0, 2] with 10001 grid points, and 



10 



x=0 





FIG. 8: The orbit along the spike point x — projected onto 
the (E+ , plane for the perturbed second (left) and third 
(right) order spike simulation. The orbits follow the false and 
the true spike orbits respectively. The initial (at T — 0) and 
final (at T — 12) points are marked by the letters i and f 
respectively. 



The initial data for a perturbed third order spike so- 
lution at r = is also given recursively: 



sN- 



X2' 



7V_ = -sS_ 



2' 



(66) 

Vx2, 

(67) 
(68) 



where (S_, A^x , Sx , N-)2 are the perturbed second order 
spike solution in above, c and s are given by 



fi 



2/2 



/l + l 



(69) 



/2 = 



Q3_ 

2Qo 



(70) 



where fi is given in (|62p . and x ^ X + Xzoom- ^2 is 
specified by numerically evaluating the constraint 

S2 = (S2)oe^«'^"/^^''". (71) 
We perforin a simulation with the parameters 



e = lO"-', 
(S2)o = 10" 



w = 9.5, 
Q2 =0, 



Q3 = o, 



(72) 
(73) 



and Xzoom — 0, X ^ [0, 2] with 10001 grid points, and 
Te[0,12]. 

The center or inner part of a perturbed second order 
spike evolves into a false first order spike, as suggested 
by Figures [H and [TOl False spikes are merely a spiky 
representation of the vacuum Bianchi type II solution. 
The center or inner part of a perturbed third order spike 



evolves into a true first order spike, as suggested by Fig- 
ure 111! The outer parts of the perturbed spikes move 
beyond the domain of simulation and are suspected to 
evolve into first order spikes, with a moving spike point. 
A global numerical scheme might be needed to follow 
their evolution, but one with enough numerical resolution 
would take months to run, which is impractical. Figure[5] 
shows that the orbit for a perturbed second (third) order 
spike later follows the predicted orbit for the false (true) 
first order spike. This suggests that higher order spikes 
break into separate first order spikes, and therefore are 
not attractor themselves. 

The approximate values for the w and u parameters 
when near a Kasner point (except the first point) for 
the X = orbits in Figure [8] are given below (rounded 
to 4 decimal points). For the left figure, two false spike 
transitions and (Ex ) curvature transition link the Kasner 
epochs. 



w 



_6.5005 f'^'^^''^'^ 4.4161 ^"-^"'^ 2.4980 



u w 2.7502 



false spike 



1.7081 



0.5003 



1.3351 



(74) 

3.0024 (75) 



Recall that false spike transitions and curvature transi- 
tion are physically the same. 

For the right figure, two frame transitions and a spike 
transition link the Kasner epochs. 

w « -5.5004 5.4999 -1.4979 1.5002 



u« 2.2502 



2.2500 ^ 4.0169 



CONCLUSION 



(76) 

3.9984 (77) 



We have found numerical evidence (from both per- 
turbed solutions and generic initial data) that the spike 
solution is part of the generalized Mixmaster attractor. 
We have found that the second and third order spikes are 
not part of the attractor, and conjecture that all higher 
order spikes are not part of the attractor. 

We summarize the above conjectures as follows: 

1. Spike transitions are a new type of oscillation on 
approach to the singularity, with each transition 
approximated by a spike solution. A spike transi- 
tion has a map oi u ^ u — 2 and is different from 
the previously known Mixmaster oscillation, which 
has a map of u — > u — 1. It occurs in a causal 
neighborhood of special 2D surfaces of worldlines 
in generic spacetimes. 

2. Higher order spike transitions (with maps u — > m — 
3, etc) split into first order spike transitions and so 
are not general, i.e. the generic behavior towards 
singularity is either u— fu— 1 or u—^u — 2. 
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Numerical solutiori Exact solution Difference 




FIG. 9: The Weyl scalar invariants (normalized by Hubble scalar) for the numerical solution, the matching exact false spike 
solution (with w — 5.5005, tq — 2.3329), and their difference, during the first false spike recurrence over the time interval [1, 3]. 



We have used symmetric data in order to hold the spike 
point fixed, so that we can zoom in on it. We believe that 
for non-symmetric data, in which the spike point can 
move (by a little when the spike is active, and sometimes 
by a lot when the spike is smoothed out), the above con- 
clusion should also hold. This remains to be confirmed 
numerically. At present we do not know how to zoom in 
on a moving spike point. 

What remains unanswered is the following. Because 
we simulate only the neighborhood of a spike, we do not 
know what happens outside this domain. We also do not 
know what happens to new spike points that are created 
and move out of the domain, how they interact with other 
spike points or false spike points. Existing numerical sim- 
ulations from the global view suffers from expensive re- 
sources needed to resolve spikes, which severely limit the 
length of simulation. We envision a new way to simulate 
spikes, by combining the zoom-in view with the global 
view. The biggest benefit of such a combination is much 
longer simulations. The zoom-in view can also provide 



boundary conditions, so that the assumption of spatial 
periodicity can be dropped. Implementing the combina- 
tion will be challenging. 

APPENDIX A: MATCHING WITH EXPLICIT 
SOLUTIONS 

For the purpose of matching the numerical solutions 
with explicit solutions, we will need the explicit spike so- 
lutions with generic time and space constants. To restore 
these constant, perform the transformation 

2 

T^T-To, X ^ ———{X - Xq). (Al) 
i-C'l jo 

The expression of the metric, the governing equations, 
and the solutions will change accordingly. In particular, 
Ei^ is now given by 

E,^ = {E.^e^-^o. (A2) 
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Numerical solution Exact solution Difference times 1000 




FIG. 10: The Weyl scalar invariants (normalized by Hubble scalar) for the numerical solution, the matching exact false spike 
solution (with w = 1.5006, to — 7.6904), and their difference, during the second false spike recurrence over the time interval 
[6,10]. 



The spike solution is now given by 



P = 2(t - To) + ln[sech(u;(T - To))] - Hf + 1] 



- ln(2Qo) 

Q = -(3ow[2(wtanh(w(r - To)) - {x - xq)) 

A = -4 ln[sech(u;(T - tq))] + 2 \n[f + 1] 

- (w^ +4)(T-ro) + A2, 



(A3) 

2 



(A4) 
(A5) 



where 



f = we'^ ^°sech(w(r - To)) [x - xq) (A6) 

Jo 



is the the factor Qe^ for the vacuum Bianchi type II so- 
lution. Correspondingly, the /3-nornialized variables be- 



come 



^ (1 + T^qrY tanh(^«(T - To)) - 1]) (A7) 



2/ w 



p + 


1 V3 


p- 


1 w 


p + 


1 V3 


V 


1 


p + 


1 V3 



sech(w(T - To)) (A8) 
sech(w(T - To)) (A9) 
[1 - wtanh(u;(T-To))]. (AlO) 



The other solutions are similarly restored. 

We take this opportunity to correct errors in the 
third minus sign in Equation (28) should be a plus sign, 
and the factor 4 in Equation (34) should not be there. 

In order to match with an explicit spike solution, we 
will need to guess the value of the parameter w. This 
can be done in two ways. The first way is to choose 
a predetermined value, the second is to obtain a guess 
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FIG. 11: The Weyl scalar invariants (normalized by Hubble scalar) for the numerical solution, the matching exact spike solution 
(with w — —3.5002, to = 4.4985), and their difference, during the first spike recurrence over the time interval [3, 6]. 



from the numerical solution. To do so we compute the 
expression 



arcsinh 



along X = 0, which equals 




w{t - To) 



(All) 



(A12) 



for the spike solution. The value tq is then obtained 
through interpolation. 

We then compute the expression 



2 

V3 



(A13) 



along X — 0. For the spike solution this expression 
equals /3. In practise the numerical solution will give 
an close-to-constant time function of w, from which we 
choose one value. For example we can take the maximum 
value of this time function. 



APPENDIX B: OBTAINING THE BKL 
PARAMETER u FOR THE KASNER EPOCHS 



Matching the Kasner epochs with Kasner solutions is 
straightforward. Recall from equation (18) of Q that for 
a Kasner solution 



w 



V3' 



(Bl) 



where w is a constant. One then obtains the local maxi- 
mum and minimum values for E_ along a worldline and 
convert them to w. Then one computes the BKL param- 
eter u from w using the following formula. 



2 
2 

l-\w\ 



\w\>3, 
1< |w| < 3, 
< w < 1. 



(B2) 
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APPENDIX C: THE WEYL SCALAR 
INVARIANTS 



where N+ = VSiV-, r = -3(iVxS- - A^-^x)- The four 
Weyl scalar invariants are computed as foUows: 



The orthonormal frame components Cabcd of the Weyl 
tensor can be conveniently expressed in terms of the elec- 
tric and magnetic components Eap and Ha^ [l5| : 

CaO/30 = Eaf3, Ca/SjS ~ —^'^ afl^jS-^t^i^ ' (CI) 

CafijO — ^^apH-^^i, (C2) 

which are then normalized by 2>(3^: 

and further decomposed as follows: 

(-2E+ V3£3 \ 
fa/3 = f + + ySf- V35x (C4) 

Vv^f2 %/3fx f+ - \/3£_/ 

and similarly for Tiafj- The components are given by 

8+ = iS+ - + S^x) + |(A^^ +iV2) + (C5) 
= i(l - 3S+)E_ + |iV+iV_ 

+ i(e-4^i?iiax-r)iVx+5i5S^ (C6) 
£x -i(l-3S+)Ex+fiV+iVx 

-i(e-4^i?iiax-r)iV- (C7) 

(C8) 

(C9) 

7i:+ = -iV-E- - TVxS- (CIO) 
H- = -S+A^_ - |A^+S_ - ^{e'^^Ei^dx - r)Sx 

(Cll) 

(C12) 
(C13) 
(C14) 



f3--75SxS2 

£:2 = i(i + \/3)S-)S2 



Hx = -S+7Vx - I^V^Sx + Ue^^E.^dx - r)S_ 



W3 = -;^iVxS2 

73^ 



CafccdC"'^^'^ = SiE^pE'^^ - H^pH'^P) (C15) 

Cabcd* C""^'* = IGE^f^H^P (C16) 

Cab^'^Ce/^Ce/"'' = -16{EjEfi''E^^ - 3EjHp''H^^) 

(C17) 

(C18) 

where *Cafccrf = ^Vab"^ Cefcd, and 77'^^'='' is the totally an- 
tisymmetric permutation tensor, with r/^^^^ = 1. 

The drawback of plotting the Weyl scalars for spikes 
is that the blow-up of the Weyl scalars towards the sin- 
gularity makes the spiky structures invisible. For ex- 
ample, in Figure 8 of Q, level curves have to be plot- 
ted to make the structure visible. In this paper we plot 
Hubble-normalized Weyl scalars so that spiky structures 
are clearly visible. The Weyl scalars are normalized as 
follows: 



/^abcd 
(3i?2)2 ' 

cd/^ ef ab 
'^ab ^cd ' L'e/ 

(3fl^)3 ' 



^abcd ^ 

(3i72)2 ' 

^ab '^cd ' t-'e/ 



2^3 



(C19) 
(C20) 



n2 = ^^N^^2 
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